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Abstract .  The  paper  addresses  the  question  of  the  optimal  selection  of  the 
shape  functions  for  p-type  finite  elements  and  discusses  the  effectivity  of 
the  conjugate  gradient  and  multilevel  iteration  method  for  solving  the 
corresponding  linear  system. 


1 .  Introduct Ion. 


A  natural  way  to  treat  the  finite  element  method  on  a  parallel  computer 
is  to  associate  to  every  available  processor  a  superelement.  The  main  compu¬ 
tational  work  is  then  on  the  superelement  level  and  car  be  made  completely  in 
parallel.  The  global  part  of  the  computation  which  is  not  well  paral lei izable 
is  then  a  relatively  minor  part  of  the  entire  computation. 

The  superelement  can  be  for  example  constructed  by  the  p-version  of  the 
finite  element  method,  by  the  domain  decomposition  technique  where  the  super¬ 
element  consists  of  a  grid  of  finer  elements,  etc. 

The  complexity  of  such  a  superelement  approach  (based  on  the  p-version) 
depending  on  the  computer  architecture  was  studied  in  [1],  [2],  See  also  [3], 
[4). 

We  can  understand  the  superelement  as  a  linear  space  S  of  trial  func¬ 
tion  spanned  by  certain  basis  (shape)  functions.  The  convergence  and  accuracy 
of  the  finite  element  method  then  depend  on  the  approximation  properties  of  S 
and  is  of  course  independent  of  the  particular  choice  of  the  basis  of  S.  On 
the  other  hand,  the  computational  effectivity  of  the  finite  element  method 
depends  also  for  example  on  the  effectivity  of  an  iterative  solver,  the  sensi¬ 
tivity  of  a  direct  solver  round-off  errors,  etc.,  which  depends  directly  on 
the  choice  of  the  basis  functions. 

The  main  mathematical  tool  of  analyzing  the  performance  of  the  numerical 
method  is  the  asymptotic  analysis  with  respect  to  some  parameters,  such  as 
p— in  the  p-version,  or  h — >0  in  the  h-version.  Such  an  analysis  is 
effective  in  practice  only  if  the  asymptotic  range  encompasses  the  practical 
values  of  these  parameters.  Unfortunately,  this  is  not  always  the  case.  Let 
us  mention  for  example  the  p-version.  Here  the  performance  of  certain  itera¬ 
tive  method  (see  below)  can  be  characterized  by  a  parameter  H(p)  which  is 
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known  to  grow  asymptotically  as  K(p)  ~  (lgp)  •  However,  as  in  practice  one 
mostly  has  p  <  10,  the  asymptotic  analysis  neglecting  the  precise  behavior 
of  #(p)  for  small  p  is  not  sufficient.  Although  the  asymptotic  analysis 
is  still  important,  practical  insight  in  the  performance  of  a  method  can  be 
obtained  only  together  with  a  proper  computational  analysis.  The  aim  of  this 
paper  is  to  make  such  a  computational  analysis  and  to  get  an  insight  into 
various  aspects  related  to  the  choice  of  basis  (shape)  for  p-type  finite  ele¬ 
ment  functions.  We  will  restrict  ourselves  here  only  to  the  case  of  a  scalar 
equation  in  two  dimensions  and  to  a  quadrilateral  element. 

A  more  general  and  thorough  analysis  will  be  given  in  forthcoming  papers. 
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2.  The  quadrilateral  superelement . 

We  will  discuss  the  case  of  a  unit-square  superelement,  which  c an  also  be 
understood  as  a  master  element. 

Let  us  consider  the  unit  square  Q  (the  master  element)  as  shown  in 
Figure  2. 1 


Figure  2.1.  The  master  element 

with  the  vertices  and  sides  r^.  Denote  by  the  set  of  shape  func¬ 

tions  used  in  the  finite  element  method.  Then  these  functions  are  of  three 
groups  associated  with  the  vertices,  sides  and  interior  of  the  element.  They 
are  called 

1)  nodal  shape  functions 

2)  side  shape  functions 

3)  internal  shape  functions. 

1)  The  nodal  shape  function  is  associated  to  a  vertex  of  the  ele¬ 

ment  n.  It  is  zero  on  the  opposite  sides  of  the  vertex  it  is  associated  to. 

2)  The  side  shape  function  is  associated  to  a  side  and  is  zero  on 

all  three  other  sides  of  the  element. 

3)  The  internal  shape  function  is  zero  on  all  four  sides,  i.e.,  it  has  a 
character  of  the  "bubble"  function. 
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In  addition  to  the  above  mentioned  properties  of  the  shape  functions  we 
require  typically  that  the  basis  is  invariant  under  rotations  of  the  coordi¬ 
nates  by  kx90°,  k  =  1,2,3.  We  require,  of  course,  that  the  span  S  of  the 
shape  function  have  good  approximation  properties,  too. 

The  shape  functions  are  not  uniquely  determined  by  the  above  require¬ 
ments.  For  example,  we  can  add  to  any  nodal  or  side  shape  function  an  inter¬ 
nal  shape  function  and  still  preserve  all  the  properties  mentioned  above. 

Usually  the  set  of  shape  functions  is  a  one  parameter  set  (say  with  the 
parameter  t  =  1,2,...).  In  the  case  of  the  p-version  t  =  p  and  in  the 
h-version  we  take  t  =  We  can  now  impose  an  additional  constraint  on  the 
set  of  shape  functions,  namely  the  requirement  that  the  shape  function  are 
hierarchic.  By  this  we  mean  that  if  is  the  basis  of  S(t)  for  given 

t  then  +  D  {1/1^},  l.e.,  the  set  +  of  shape  functions  can 

be  obtained  by  an  augmentation  of  the  set  The  property  of  hierarch- 

icity  is  desirable  if  the  computation  is  to  be  made  with  many  values  of  the 
parameter  t,  such  as  in  adaptive  computations. 
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3.  Examples  of  shape  functions. 

In  this  section  we  will  introduce  some  sets  of  the  shape  functions  that 
we  will  consider  in  the  next  sections. 

A.  The  hierarchic  set  of  the  p-verslon.  This  set  is  for  example  used 

in  the  computer  program  PROBE. 

1)  The  nodal  shape  functions.  There  are  the  usual  bilinear-  functions: 

fyc.T))  =  I(i-oa-n) 

N2(C,t))  =  i(l+<;)(l-7j) 

n3  ( c .  -n }  »  j(l+$)(i+n) 

N4(€.7j)  =  i(l-?)(l+T)) 


2)  The  side  shape  functions.  There  are  ( p— 1 )  shape  functions  associ¬ 
ated  with  every  side  v  =  1,2, 3, 4.  These  are  defined  as 

n!1](£,t})  =  1(  l-r,)*i  (^) .  i  =  1,2 . p-1, 

f  21  1 

Nj  (?,T))  =  i(l+?)*i(T»),  i  =  1,2 . p-1, 

n[31(S,t,)  =  i  =  1,2 . p-1, 

n[41(?,t))  =  i  =  1,2 . p-1. 

where 


*!<*) 


i-1 

2 


A 

P^tldt 

-1 


and  P  j  ( t ) 
needed  in 


is  the  Legendre  polynomial  of  degree  j.  The  term  (-1)*  is 
N^3^  and  to  obtain  invariance  with  respect  to  the  rotation 
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of  coordinates. 


3)  The  internal  shape  functions.  For  p  <  4  there  are  no  internal 
shape  functions.  .For  p  >  4  there  are  (p-2)(p-3)/2  internal  shape  func¬ 
tions  defined  as 

N^(C.T))  =  (l-^Hi-i^lP^OPjfn),  0  <  i  +  j  <  p-4. 

For  example,  if  p  =  8  there  are  47  shape  functions  consisting  of  4  nodal,  28 
side  and  15  internal  shape  functions. 

The  set  Q'  is  the  minimal  set  which  Includes  all  polynomials  of  degree 
P 

p  and  which  has  the  properties  listed  above.  This  guarantees  the  good 
approximation  properties  of  the  span.  Let  us  mention  that  this  set  does  not 
include  all  polynomials  of  degree  p  in  both  variables  £  and  rj  separate¬ 
ly.  The  span  is  identical  with  the  span  of  the  serendipity  elements  and  is 
denoted  by  as  in  [5].  The  introduced  set  of  the  shape  functions  is 

obviously  heirarchal  one  with  respect  to  the  parameter  p  =  degree.  These 
shape  functions  are  sufficiently  orthogonal  in  typical  applications,  so  that 
their  use  does  not  lead  to  large  problems  with  round-off  errors  in  connection 
with  direct  solvers. 


B.  The  hlerarchal  set  of  the  p-version.  This  set  includes  all  polyno¬ 

mials  of  degree  p  in  each  variable.  Hence  it  differs  from  the  set  Q'  only 

P 

2  2 
by  using  now  (p-1)  Internal  shape  functions.  The  set  Qp  has  thus  (p+1) 

2 

shape  functions  in  comparison  with  -- + P  +  3  shape  functions  in  the  set 


C.  The  quasi  orthogonal  set  0^  of  the  p-version.  Let 

i  =  l,...,p-l,  be  polynomials  of  degree  <  p,  such  that  4>i(-l)  =  4^(1)  = 
0,  and 
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+  1 

*'($)*j(€>d$  =  5lJ.  i.  j  =  1 . P-1 

-1 

+  1 

(C)*jCO<lC  =  0.  i  *  J. 

J_1 


This  set  is  uniquely  determined  as  the  (finite  element)  eigenfunctions  of  the 
eigenvalue  problem:  Find  the  pairs  (4n,A^)  e  V^xR,  i  =  l,...,p-l,  such 
that 


,+  l 


<t> '  ( C ) v '  (£)d£  =  A 

-1 


*  vd$ 


-1 


V  v  €  V 

P 


where 


V 

P 


{v|v  is  a  polynomial  of  degree  <  p  on  [-1.1]  and  v ( —  1 ) 


v ( 1 )  =  0}. 


The  eigenvalues  A^  are  positive  and  distinct.  Further,  let  ^(£),  i  = 

l,...,p-l,  be  another  set  of  polynomials  of  degree  <  p  such  that  </>.(- 1 ) 

=  1 ,  <p .  ( 1 )  =  0  and 
l 


+  1 


,1 


?;v'd<;  +  \l 


-1 


'fit vd^ 


=  0  V  V  e  V 


where  A^  are  the  eigenvalues  mentioned  above.  The  functions  <p ^  are 
uniquely  determined.  To  see  it  we  w’-lte 


?*(<;)  =  X(€)  +<fi\(S). 


x<€)  =  gd-O- 


Then  <pA$)  €  V 
i  P 


is  the  (finite  element)  solution  of  the  problem 
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,1 


v'd?  +  Al 


+  1 


d?  =  - 


( v'  +X,  *v)d£,  v  e  V  , 
1  P 


and  since  >  0  y  exists  and  is  uniquely  determined.  Now  we  define 

1)  the  nodal  shape  functions.  They  are  the  same  as  in  the  set  Q' . 

P 

2)  The  side  shape  functions.  The  side  shape  functions  associated  with 
T  are 

n[  1 1  (C.  -o)  =  *i(Ovi(u).  i  =  1.2 . p-1. 

The  shape  functions  associated  with  the  sides  T.,  j  =  2,3,4,  are  defined 

3 

analogously.  Hence  as  before  we  have  (p-1)  side  shape  functions  associated 

with  every  side  T  .. 

J 

2 

3)  The  internal  shape  functions.  There  are  (p-1)  internal  shape 
functions  as  in  the  system  Q^.  The  shape  functions  associated  with  T ^  are 
then 

Nk°i(S’v)  =  *k(?)V7,)'  k,(- =  1 . p_i- 

This  system  0p  has  all  properties  as  the  system  except  that  the  system 

0^  is  not  a  hierarchal  one.  On  of  the  hand,  it  has  various  useful  orthogo¬ 
nality  properties  which  will  be  mentioned  later. 


D.  The  h-version  set  of  the  shape  funct ions.  Here  the  span  S  will  be 

2 

the  set  if  piecewise  bilinear  functions  on  the  uniform  mesh  with  h  =  -.  The 

P 

shape  functions  are  the  usual  "hat"  functions  which  obviously  car  be  divided 

into  the  groups  of  the  nodal  side  and  internal  shape  functions. 

In  addition  we  will  consider  the  set  H*  where  the  nodal  shape  functions 

of  the  set  H  are  replaced  by  those  of  the  set  Q'  . 

P  P 

E.  The  trigonometric  set  T  of  the  shape  functions.  For  theoretical 

p - 
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reasons  we  will  also  consider  the  set  of  shape  functions  where  only  side  shape 
functions  are  present  and  defined  by 

n[13(£,t))  =  cos  $  sin  h  i  =  1,2 . p-  1, 

for  side  T  and  analogously  for  the  other  sides.  Note  that  these  functions 
are  harmonic. 
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4.  The  optimal  selection  of  the  shape  functions. 


Let  the  set  {< of  the  shape  functions  be  given.  Then  in  the  finite 

element  method  we  construct  first  the  local  stiffness  matrix  A  =  {a.  .} .  a.  .  = 

- xj  1J 

where  B(u,v)  is  the  bilinear  form  on  which  the  finite  element 
method  is  based.  The  global  stiffness  matrix  is  constructed  by  the  usual 
assembly  process.  In  this  paper  we  will  restrict  our  analysis  to  the  case  of 
one  superelement  only.  We  will  assume  that  the  bilinear  form  is  associated  to 
the  Laplace  operator,  i.e., 


(4.  1) 


B (0.  ,0.) 

i  J 


'30i 

dT 


n<- 


30i 

3tT 


J 


d£d7}. 


The  stiffness  matrix  A  is  uniquely  (up  to  permutation)  defined  by  the  set  of 
shape  functions.  We  will  always  assume  that  the  shape  functions  are  rescaled 
(preconditioned)  so  that  a  =  1  (i.e.,  all  diagonal  terms  of  A  are  equt  1 

*  »  A 

to  one).  This  will  sometimes  be  called  trivial  preconditioning. 

Now  we  can  formulate  what  we  mean  by  optimal  shape  function  selection 
according  to  various  criteria. 


A.  The  criterion  of  the  minimal  condition  number. 

We  will  say  that  the  shape  functions  are  optimal  if  the  condition  number 
of  the  local  stiffness  matrix  (after  trivial  preconditioning)  is  minimal  among 
all  choices  of  the  shape  functions  (with  the  same  span)  which  preserve  the 
categories  of  the  nodal,  side  and  internal  shape  functions,  the  invariances 
with  respect  to  the  rotation  defined  earlier.  We  can  also  impose  additional 
constraints  as  the  hierarchy  of  of  the  shape  functions. 

The  motivation  of  this  notion  is  to  get  the  most  effective  conjugate  gra¬ 
dient  method  (see  section  6)  or  to  maximize  the  numerical  stability  of  the 
direct  solution  method.  The  exact  structure  of  these  optimal  shape  functions 


10 


is  not  known. 


Realizing  that  the  performance  of  the  conjugate  gradient  method  depends 
not  only  on  the  condition  number  (e.g. ,  on  the  clustering  of  the  eigenvalues) 
we  have 

B.  The  criterion  of  the  maximal  decay  number  (i.e. ,  reduction  of  the  residuum 
by  one  iteration)  in  the  conjugate  gradient  method. 

Here  the  obvious  goal  is  to  achieve  the  required  accuracy  with  minimal 
number  of  iterations. 

We  can,  of  course,  design  other  criteria  of  optimality.  In  this  paper  we 
will  mostly  address  the  criterion  A.  We  approach  the  optimization  problem 
formulated  above  via  numerical  experiments  so  as  to  obtain  the  information  in 
the  practical  range  of  p,  say  for  1  <  p  <  20. 
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p 


We  will  first  analyze  the  condition  number  of  the  (local)  stiffness 

matrix  A  when  using  the  bilinear  form  (4.1).  Since  the  first  eigenvalue 

A 

A,  =  0  we  define  the  condition  number  Jf(A)  =  -r —  and  call  A  and 
1  Ag  max 

A  .  =  A„  the  dominant  eigenvalues.  To  get  a  proper  insight  we  will  compute 

min  2 

the  condition  number  of  various  combinations  of  the  shape  functions.  The 
results  are  shown  in  Table  5. 1. 


Table  5. 1. 

Condition  number  H  of  various  portions  of 

the  stiffness  matrix  of  set  set  Q' . 

P 


p 

Internals 

+ 

Sides 

Nodal s 

Internals 

Sides 

Sides 

+ 

Nodals 

I:  .  ~nals 

Nodals 

Internals 

+ 

Sides 

1 

1.50 

1.50 

1.50 

2 

5.  78 

1. 36 

5.77 

1. 50 

1. 36 

3 

5.78 

2.33 

5.77 

1.50 

2.33 

4 

22.  1 

1.00 

2.70 

5.96 

1. 50 

22.  1 

5 

23.9 

1.00 

2.84 

5.96 

1.50 

23.6 

6 

62.8 

3.09 

2.92 

5.  96 

3.09 

62.8 

7 

67.7 

4.53 

2.95 

5.96 

4.53 

66.7 

8 

203.0 

12.8 

2.97 

5. 96 

12.8 

203.0 

The  first  column  in  Table  5.1  depicts  the  condition  number  of  the  whole 

3 

stiffness  matrix.  The  condition  number  H(A)  grows  as  p  (see  Figure  5.1). 
Upon  comparing  the  columns  of  Table  5.1  we  see  that  the  dominating  factor  in 
the  growth  of  the  condition  number  is  the  coupling  between  the  side  and 
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Figure  5.1.  The  condition  number  K(A)  for  the  set  Q^. 

internal  shape  functions.  Other  coupling  such  as  sides-sides,  sides-nodals, 
internals-internals,  are  relatively  much  weaker  except  in  the  low  range  (1  < 
p  <  3)  where  there  are  no  internals.  We  are  then  led  to  the  conclusion  that 
in  seeking  for  an  optimal  set  of  shape  functions,  the  first  step  should  be  to 
reduce  this  coupling  between  the  side  and  internal  shape  functions.  Below  we 
will  remove  the  coupling  completely  by  eliminating  (condensation)  the  inter¬ 
nals  (this  can  be  done  in  parallel  when  many  superelements  are  present). 
Although  this  very  likely  does  not  lead  to  an  exactly  optimal  basis,  we 
conjecture  that  we  can  in  this  way  construct  shape  functions  which  are  close 
to  the  optimal  ones  (in  the  sense  of  criterion  A). 

The  above  mentioned  elimination  can  be  interpreted  as  changing  the  side 
shape  function  so  that  they  are  "harmonic"  in  the  finite  element  sense.  This 
removes  the  coupling  between  internals  and  side  shape  functions  but  increases 
the  coupling  between  side  shape  functions.  Of  course,  the  hierarchy  of  the 
shape  functions  is  lost.  Note  that  the  nodal  shape  functions  are  bilinear  and 
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hence  harmonic,  so  they  are  orthogonal  to  the  internal  shape  functions. 

We  will  now  address  this  decoupling  in  more  detail.  Given  p  we  will 
order  the  shape  functions  so  that  the  internals  come  first,  i.e.,  we  write 


¥  =  [#tll,*[2I]T 


where 


r .  [  1 1  .m, 

*  =  . 1 


are  the  internal  shape  functions.  The  stiffness  matrix  can  be  written  now  in 
a  block  way 


A  = 


A11  A12 

*L>  *22 


Here  the  block  corresponds  to  the  internals.  The  decoupling  is  per¬ 

formed  by  a  transformation  into  the  new  basis 

#  =  ,¥[21 ] 

where 


(5.  1) 


i121 


Here  M 2  is  a  diagonal  (normalizing)  matrix  to  be  defined  shortly.  In  this 
new  set  of  shape  functions  the  stiffness  matrix  takes  the  form 


(5.2) 


where 


(5.3) 


A  =  BAB  = 


11 


B  = 


fl,  o  1 

1 

0  M_ 

L  2J 

0  "2*22*2 


1 

«  X2 
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is  the  Schur  complement 


and 


(5.4) 


T  -1 

A22  =  A22  "  A12A11A12' 


By  1^  and  we  denote  the  identity  blocks,  k^ 


of  A  with  respect  to  A 


11' 


We  finally  define  the  diagonal  matrix 


M  = 


ri  ° 


so  that  [A]  ,  =  1  for  all  i,  i.e.,  we  set 

A  »  1 


(5.5) 


lM,u  * 


Since  the  nodal  shape  functions  are  bilinear  (and  harmonic)  there  is  no 

coupling  between  the  nodals  and  interiors.  Therefore  the  transformation  (5.1) 

*[2] 

changes  only  the  side  shape  functions.  The  shape  functions  i/i.  are  such 
that 

-  f  2  1 

B  (*[,*)  =  0 

for  any  <p  =  0^ .  i.e.,  are  discretely  harmonic. 

If  AjjAi2  in  (5.4)  is  computed  using  the  Cholesky  decomposition  A^  = 

T  ~  [  l  ] 

LL  then  we  obtain  as  a  by-product  new  interior  basis  functions  <p  which 

" [ 1 ]  -Til] 

are  orthogonal,  i.e.,  with  0^,  =  L  0^  we  get 


B(i[1].ij1])  =  0[1]tl_1al  =  5itJ, 


and  hence  the  stiffness  matrix  for  these  new  internals  is  the  unit  matrix. 

-  [ 1 1  * [ 2] 

The  system  of  the  shape  functions  0  =  [0  ,  0  ]  constructed  above  will  be 

*  * [ i ]  - [2] 

denoted  by  Q'  ,  and  the  system  of  shape  functions  {0  ,  i l>  }  will  be 

P 

_  • 

denoted  by  Q' . 

P 

Remark  5 . 1 . .  Note  that  k^  is  the  stiffness  matrix  that  remains  in  the 
finite  element  system  where  the  degrees  of  freedom  for  the  internal  shape 
functions  are  eliminated  (condensed  out).  We  will  call  the  element  where  the 
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interior  degrees  of  freedom  are  removed  a  condensed  element.  Thus  in  the 


condensed  element  there  are  no  internal  shape  functions.  The  side  shape  func¬ 
tions  are  "discretely  harmonic"  with  respect  to  the  original  internals,  and 
the  condensed  matrix  is  the  Schur  complement  of  the  original  stiffness  matrix 
with  respect  to  its  interior  block. 

Let  us  now  consider  matrix  A  as  defined  by  (5.2)  through  (5.5).  In 

Table  5.2  we  give  the  dominant  (i.e.,  the  smallest  non-zero  and  the  largest 

- 

eigenvalue  of  and  A^2  =  ^^22^2 ^  33  we**  33  t*ie  cond^tion  number  A 

for  p  =  1, ....  8. 


Table  5.2. 


Dominant  eigenvalues  of  A^  and  -  ^2^26  and  H(A) 


for  the  shape  functions  of  Q' 

P 


p 

A  ,  (A,,) 
min  11 

\ln(A22> 

A  (A..) 

max  11 

A  (A__) 

max  22 

K(A) 

1 

1.0 

1.50 

1 . 50 

2 

0.346 

2.00 

5.78 

3 

0.346 

2.00 

5.78 

4 

1.0 

0.  253 

1.0 

2.20 

8.69 

5 

1.0 

0.242 

1.0 

2.24 

9.26 

6 

0.524 

0.  190 

1.62 

2.26 

11.9 

7 

0.402 

0.  177 

1.82 

2.26 

12.8 

8 

0.  183 

0.  141 

2.35 

2.  26 

16.7 

Comparing  Tables  5.1  and  5.2  we  see  that  removing  the  coupling  between 
the  internal  and  the  side  shape  functions  reduces  the  coiwuion  numoei  of  the 
stiffness  matrix  quite  dramatically.  We  see  that  H{A)  grows  with  p  much 
more  slowly  than  H(A).  Note  also  that  since  A^  has  no  contribution  to 
K(A)  in  view  of  Table  5.2,  there  is  no  reason  for  constructing  better  (e.g.. 
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orthogonal)  internal  shape  functions  \fj  at  least  in  the  range  1  <  p  <  8. 

•  _  • 

In  this  range  the  condition  numbers  of  the  set  Q'  and  Q'  are  the  same 

P  P 

and  also  the  condition  number  of  the  condensed  element  is  the  same  as  that  of 


So  far  we  have  not  addressed  the  shape  functions  on  the  sides  r  .  To 

assess  their  influence  on  the  condition  number  let  us  construct  orthogonal 

'  [2] 

normalized  shape  functions  ^  and  the  corresponding  stiffness  matrix  A 

A 

A 

namely  such  that  fA]^  =  0  for  *  *  j  *n  following  cases: 

-[2]  ~[2] 

a)  \j).  and  are  both  internal  shape  functions; 

^  J 

"[2]  "[2] 

b)  is  an  internal  shape  function  and  is  the  side  shape 

funct ion; 


A[2]  A[2] 

c)  i /n  and  \(i ^  are  both  side  shape  functions  corresponding  to  the 

one  side. 


These  shape  functions  are  obviously  decoupled  in  a  maximaly  orthogonal 
way  under  the  constraint  imposed  in  section  2,  namely  preserving  the  groups  of 
nodal  side  and  internal  shape  functions. 

We  can  construct  the  shape  functions  els  follows.  Let  ^  be  the  shape 
functions  already  constructed.  Let  be  the  matrix  (of  the  same  size  as 

the  matrix  A)  such  that 


1) 

tA0]ii  = 

1 

for  all  i; 

2) 

II 

H-t 

<  < 

a 

(A] 

^  for  i  *  J  in  any  of  the  cases  a,b,c  mentioned 

above; 

3) 

[Vij  = 

0 

otherwise. 

Let 

further 

O 
<  < 

T 

=  CC  be  the  Cholesky  decomposition  of  Aq.  Then  we 

define  the  new  shape  functions 
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c 


The  corresponding  stiffness  matrix  is  then 


(5.6) 


~  -I-  -T 

A  -  C  AC  . 


These  shape  functions 
stiffness  matrix  A  (with 

•  m 


be  denoted  by 


% 


0  satisfy  all  desired  properties  as  well  as  the 

[A].  .  =  1).  This  set  \p  of  shape  functions  will 
1  1 


Remark  5.2.  For  symmetry  reasons  it  is  obviously  sufficient  to  compute  the 
new  shape  functions  associated  to  only  one  side,  say  r.,  and  then  to  use 
simple  transformations  to  get  the  side  shape  functions  on  the  other  sides. 


Remark  5.3.  The  above  method  of  constructing  maximally  orthogonal  functions 

is  universal  in  the  sense  that  it  obviously  apply  to  any  set  of  shape  func- 

**  «* 

tions.  For  example,  we  can  construct  in  a  similar  way  the  sets  Q  ,  H  and 

P  P 

+•« 

11^  We  can  construct  the  shape  functions.  There  are  many  other  ways  of 

constructing  maximally  orthogonal  systems  of  shape  functions.  For  example, 

the  set  0  mentioned  above  is  such  as  is  easily  checked. 

P 

Remark  5.4.  The  shape  functions  \p  and  i/»  are  obviously  not  hierarchic. 

For  computational  reasons  it  may  be  desirable  to  keep  the  new  stiffness  matrix 
in  the  product  form  (5.2),  (5.6),  so  that  the  non-hierarchic  part  (matrices  B 
and  C  * )  is  separated  from  the  hierarchic  part  (matrix  A).  In  the  many- 
element  case  (such  as  on  grid  of  quadrilaterals),  it  is  sufficient  to  compute 
B  and  C  1  on  a  reference  element  only.  This  sacrifices  in  general  the 
orthogonality  but  essentially  preserves  the  good  conditioning  if  the  elements 
are  not  too  much  distorted. 

In  Table  5.3  below  we  give  the  dominant  eigenvalues  and  the  condition 
number  H(A)  for  p  =  1,...,8. 
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Table  5.3. 


Dominant  eigenvalues  and  condition  number 


for  the  shape  functions  of 


% 


p 

A  ,  (A) 
min 

A  (A) 

max 

H(  A 

1 

1.0 

1.50 

1.50 

2 

0.  346 

2. 00 

5.78 

3 

0.346 

2.00 

5.78 

4 

0.  253 

2.20 

8.69 

5 

0.244 

2.  24 

9. 19 

6 

0.  200 

2.26 

11.3 

7 

0.  187 

2. 27 

12.1 

8 

0.  157 

2.28 

14.5 

Comparing  Tables  5.2  and  5.3  we  see  that  the  orthogonal izat ion  of  the 
side  shape  functions  (i.e.,  the  maximal  orthogonal ization)  does  reduce  the 
condition  number  but  the  effect  is  not  large  in  the  range  1  <  p  <  8,  (at 
most  13%). 


19 


6.  Comparison  of  various  sets  of  shape  functions . 

We  will  compare  in  this  section  the  condition  numbers  of  the  stiffness 

matrices  for  the  sets  of  shape  functions  which  were  introduced  in  the  previous 

sections.  Based  on  the  analysis  we  made  in  section  5  for  the  system  we 

restrict  ourselves  only  to  the  case  of  condensed  elements  (see  Remark  5.1). 

As  we  have  seen,  this  restriction  is  equivalent  to  the  case  where  the  sets 
~ 

Q'  ,  Q  ,  H  ,  etc. ,  instead  of  Q' ,  Q  ,  H  ,  etc. ,  are  used.  The  shape  func- 
P  P  P  P  P  P 

tion  of  the  set  T  are  harmonic  functions.  Hence  we  consider  them  as  con- 
P 

densed  shape  functions  (l.e.,  T  =  T  ). 

P  P 

In  Table  6.1  we  present  condition  numbers  of  the  condensed  stiffness 

matrices  for  various  sets  of  shape  functions  for  p  =  4n +  1 ,  n  =  1 . 5  (as 

before,  the  shape  functions  are  normalised  so  that  the  stiffness  matrix  has 
ones  on  the  diagonal). 


Table  6.1. 

Condition  number  of  the  stiffness  matrix  for  p  =  4n+l,  n  =  1 . 5 

for  various  sets  of  shape  functions. 


1 

2 

3 

4 

5 

6 

7 

p 

^  • 

Q' 

** 

q: 

* 

Q 

0 

H 

H+* 

* 

T 

P 

p 

P 

P 

P 

P 

P 

5 

9.3 

9.2 

15.5 

14.2 

5.  4 

14.7 

11.4 

9 

16.3 

14.8 

28.0 

22.  4 

10.  2 

20.  2 

15.8 

13 

23.  4 

19.5 

39.8 

28.5 

14.  9 

45.2 

18.7 

17 

30.7 

23.8 

51.3 

33.5 

19.7 

62.3 

20.9 

21 

37.8 

27.  4 

62.5 

37.8 

24.4  | 

80.4 

22.7 

We  see  that  the  condition  number  of  the  condensed  stiffness  matrix  grows 

-+• 

relatively  slowly  with  p  in  all  cases  except  for  H  It  car  be  theoreti- 

•  • 


cally  shown  that  in  the  caLses  0  ,  Q 

P  P 

2 

at  worst  of  order  O((log  p)  ),  i.e.. 


and  Tp  the  growth  is  asymptotically 
there  is  a  constart  independent  of  p 
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such  that 


(6. 1)  H(p)  <  Ct (log  p)2]. 

We  refer  to  [6],.  [7]  for  relevant  results.  Figure  6.1  shows  the  condition 

2 

number  H(p)  in  the  lx(log  p)  scale  (in  the  range  5  <  p  <  21)  for  all 
considered  ^ases.  For  example,  in  the  case  of  the  set  0^,  H(p)  can  be 
expressed  rather  accurately  by 

(6.  la)  H(p)  »  5.6  +  3. 5( log  p)2 

in  the  range  9  <  p  <  21,  or  by 

(6.1b)  K(p)  *  2. 3  +  2. 6  log  p  +  3.0(log  p)2 


in  the  who  1 e  range  5  <  p  S  21. 

_  •  _+• 

In  the  cases  Q'  ,  Q  ,  H  the  growth  of  the  condition  number  is 
P  P  P 

expected  to  be  larger  due  to  the  nonorthogonality  of  the  side  shpae  functions. 


This  nonorthogonality  is  largest  in  the  case  H 


Here  one  can  show  that  the 


growth  role  is  asymptotically  at  least  K(p)  ~  p  and  at  worst  H(p)  ~ 

2 

p( log  p)  .  The  impact  of  these  nonorthogonalities  is  clearly  visible  in 
Figure  6.1. 


Condition  number  H(p)  for  various  sets  of  shape  functions 


1)  =  Q'  2)  =  Q'  ,  3)  =  Q  ,  4)  =  Q 

P  P  P  P 


5)  =  H  ,  6)  =  H 

P  P 


7)  =  T 
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Let  us  finally  point  out  that  in  the  practical  range  of  p  (see  Table 


6.1  and  Figure  6.1)  the  element  of  the  type  Q'  seems  superior  to  Q  -type 

P  P 

element  (with  the  same  value  of  p)  even  if  the  side  shape  functions  of  Q 

are  orthogonal ized  (i.e.,  when  the  set  0  is  considered)  and  those  of  Q' 

P  P 

_  * 

sure  not  (i.e.,  when  we  use  the  set  ).  In  the  orthogonal ized  cases  the 

*  * 

condition  numbers  of  the  set  Q'  are  about  30%  less  than  of  the  set  0  . 

P 

*  * 

It  also  seems  from  Figure  6.  1  that  the  growth  of  H(p)  for  the  set  and 

2 

0^  is  the  same  (i.e.,  H(p)  *  0( log  p)  )  but  this  has  not  been  proven  theo¬ 
retical  ly  yet. 


Remark  6. 1.  It  seems  that  generally  the  condition  number  of  the  condensed 
polynomials  is  better  when  less  internal  shape  functions  are  present  in  the 
original  set.  For  example,  let  us  modify  the  set  Qp  so  that  it  includes  the 
internals  of  the  set  Q  (this  set  will  be  denoted  1  .  Q  ). 

q  p.  q 

Table  6.2  shows  the  condition  number  of  the  reduced  stiffness  matrix 
corresponding  to  two  adjacent  sides  and  the  node  inbetween  for  p  =  11  and 
varying  q.  We  see  that  the  condition  number  grows  monotonical ly  with  q. 

This  indicates  that  the  coupling  between  the  sides  becomes  stronger  as  q 
grows. 

Table  6.2. 

Condition  number  of  the  set  Q, ,  reduced  to 

1 1 ,  q 

two  adjacent  sides  and  the  node  in  between. 


q 

2 

3 

4 

5 

6 

7  ! 

8 

H(q) 

3.59 

5.82 

8.01 

10.  1 

12.  2 

14.2 

16.  2 

q 

9 

10 

11 

12 

13 

14 

15 

H(q) 

18.  2 

20.2 

22.  1 

23.0 

23.7 

23.9 

24.  1 

22 


The  results  shown  in  Table  6.2  are  in  constrast  with  the  theory  behind  the 
estimate  (6.1)  which  is  proven  only  if  q  >  p,  see  [7], 

Remark  6.2.  So  far  we  have  assumed  in  all  cases  that  the  nodal  shape  func¬ 
tions  are  bilinear.  The  practical  reasoning  behind  this  is  that  the  span  of 
these  shape  functions  contains  the  constant  function,  which  makes  it  possible 
to  use  the  nodal  block  of  the  stiffness  matrix  as  an  effective  precondition 
for  the  iterative  method  in  the  many-element  case.  For  more  see  [6],  [7], 

It  is  of  interest  to  see  how  the  couplings  between  the  nodal  and  side 
functions  affects  the  conditioning  of  the  condensed  stiffness  matrix.  In 

Tables  6.3a,b  we  give  some  results  indicating  the  strength  of  this  coupling. 

~  * 

In  Table  6.3a  we  show  the  condition  number  of  the  set  Q'  and  of  the  same 

P 

set  when  the  nodals  are  removed  from  the  original  set  (this  set  will  be 

denoted  by  Q'  1.  In  Table  6.3b  we  show  the  condition  number  of  Q  and  the 
p,  u  p 

same  set  when  the  nodal  Eire  condensed  out  (eliminated).  This  set  will  be 
denoted  by  Qp  QQ. 

The  columns  5  and  6  in  Table  6.1  indicate  how  the  condition  number  in  the 
h-version  is  affected  when  the  nodal  shape  functions  are  changed  from  the 
usual  "hat"  function  to  the  bilinear  one. 
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Table  6.3a. 


The  condition  number  of  the  sets  Q'  and  Q'  _. 

P  P.0 


•  • 


p 

5p 

Qp.O 

2 

5.78 

1.36 

3 

5.78 

2.  33 

4 

8.69 

6.00 

5 

9.26 

6.00 

6 

11.9 

8.32 

7 

12.8 

8.  13 

8 

16.  1 

11.2 

Table  6.3b. 

The  condition  number  of  the  sets  Q  and  Q 


P 

S 

Qp,00 

2 

8.69 

6.0 

3 

8.78 

6.0 

4 

15.5 

10.9 

5 

15.5 

10.9 

6 

22.0 

15.3 

7 

22.0 

15.  3 

8 

28.0 

19.6 

9 

28.0 

19.6 

10 

34.0 

23.8 

11 

34.0 

23.8 

We  see  that  in  the  practical  range  of  p  the  coupling  between  the  nodal 
and  the  side  shape  functions  may  affect  the  condition  number  by  a  factor  1.5 
to  3.  The  effect  is  strongest  in  the  case  of  the  h-version. 

We  underline  that  in  practice  when  multi -element  approach  is  used,  we  do 
preconditioning  by  p  =  1,  i.e.,  we  are  condensing  out  the  nodal  shape 
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functions  and  hence  the  second  columns  in  Tables  6.3b  essentially  governs  the 
effectivity  of  the  conjugate  gradient  method. 


Let  us  summarize  the  main  conclusions  we  have  seen  so  far. 
a)  The  main  source  of  the  large  condition  number  of  the  stiffness  matrix 
of  the  p-type  element  with  hierarchic  shape  functions  (sets  0^,(2^)  is  the 
strong  coupling  between  the  side  and  internal  shape  functions.  This  coupling 
can  be  removed  only  by  introducing  a  non-hierarchic  set  of  shape  functions. 


and 


The  orthogonal izat ion  of  the  side  shape  functions  (for  the  sets 
improves  the  condition  number,  the  more  the  higher  the  value  of  p. 


In  the  range  1  <  p  £  10  the  effect  of  orthogonal izat ion  is  at  most  40%. 

c)  The  nodal  functions  contribute  to  the  condition  number  by  about  30%. 
It  is  advantageous  to  remove  them  in  the  preconditioning  phase. 

d)  The  condition  number  of  the  condensed  matrix  improves  in  the  practi¬ 
cal  range  of  p  if  there  are  less  internal  shape  functions.  Hence  the  set 
Qp  is  preferable  over  the  set  Q^. 

* 

e)  The  condensed  h-type  superelement  (set  H^)  behaves  in  a  very 

similar  way  as  the  Q  -element. 

P 

f)  The  condition  number  of  the  condensed  Q^-  or  H^-type  element  grows 

2 

approximately  linearly  with  (log  p)  sis  p — >»  provided  that  the  side  shape 
functions  are  orthogonal.  In  the  range  IS  p  <  20,  deviations  from  this 
theoretical  growth  are  minor  for  all  considered  sets. 
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7.  Iterative  methods. 

In  this  section  we  will  consider  some  aspects  of  the  selection  of  shape 
functions  on  the  performance  of  two  basic  iterative  methods,  the  conjugate 
gradient  method  and  the  multi-level  method. 

A.  The  conjugate  gradient  method.  We  consider  the  case  of  only  one  element 
of  type  or  Qp  with  nonorthogonal  shape  functions.  The  multi-element 

conjugate  gradient  method  is  analyzed  in  [2],  [7],  where  also  various  aspects 
of  parallel  implementation  are  addressed.  Table  7.1  shows  the  average  conver¬ 
gence  role  (the  decay  factor)  of  10  iterations  using  the  condensed  matrix  when 
random  initial  solution  was  used.  More  precisely.  Table  7.1  reports  p^  = 

{energy  norm  of  the  error  after  10  steps/energy  norm  of  the  initial 

,1/10 

error) 


Table  7.1. 


The  average  convergence  rate  of  the 

_  *  • 

conjugate  gradient  method  for  the  condensed  sets  Q'  and  Q  . 

P  P 


p 

'io(V 

4 

.251 

.380 

5 

.328 

.  400 

6 

.398 

.503 

7 

.465 

.468 

8 

.419 

.577 

9 

.450 

.492 

10 

.442 

.580 

11 

.  470 

.525 

_  *  * 

The  non-monotone  behavior  of  p,.(Q'  )  and  p  (Q  )  shown  in  the  Table 

10  p  P  P 

_  * 

is  due  to  random  choice  of  the  initial  error.  We  see  that  the  set  O'  per- 

P 


forms  better  thatn 


as  expected  from  the  previous  study  of  the  condition 
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number  study  in  the  previous  section.  Note  that  for  p  <  3  the  conjugate 
gradient  method  converges  in  less  than  ten  steps  in  both  cases. 

In  Table  7.2  we  compare  the  average  convergence  rate  after  ten  iterations 
to  the  theoretical  worst-case  rate,  see  [8], 


where  H  is  the  condition  number  of  the  stiffness  matrix. 


Table  7.2. 

The  actual  average  convergence  rate  and  the  theoretical  (worst)  rate 

after  10  iterations  for  the  condensed  set  Q'  and  Q  . 

P  P 


p 

P10(5p  > 

~pior%  > 

P10(%) 

•W 

6 

.398 

.590 

.503 

.695 

11 

.470 

.679 

.525 

.758 

We  see  that  the  conjugate  gradient  method  converges  significantly  faster  than 
it  would  be  in  the  theoretical  worst  case.  It  is  likely  related  to  the  clus¬ 
tering  of  the  eigenvalues  of  the  matrix.  In  Figure  7. la,  respectively  7. lb, 
we  have  depicted  the  density  (histogram)  of  the  eigenvalues  of  the  condensed 
element  of  the  set  Qp,  respectively  the  set  0^  for  p  *  21.  We  see  that 
the  spectrum  is  sparse  at  its  ends,  most  of  the  eigenvalues  being  clustered 
around  the  point  A  =  1.0.  (For  example,  56  of  the  83  are  non-zero  eigen- 
eigenvalues  of  the  system  021  are  located  on  the  interval  [0.9,  1.1].)  The 
clustering  is  more  pronounced  in  the  case  0p  (where  the  side  shape  functions 
are  orthogonal lzed)  in  comparison  with  the  case  of  the  set  Qp.  Apparently, 
this  explains  why  the  actual  rate  of  convergence  reported  in  Tables  7. 1  and 
7.2  is  better  than  that  for  the  worst  possible  spectrum. 
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I 

Figure  7. la 


The  density  of  the  eigenvalues  of  the  set 


57 


Figure  7. lb 

The  density  of  the  eigenvalues  of  the  set  0^ 
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Let  us  finally  mention  that  the  above  qualitative  nature  of  the  spectrum, 
i.e.,  the  clustering  around  the  point  A  =  1  with  sparse  spectrum  elsewhere 
is  characteristic  to  all  elements  (i.e.,  sets  of  shape  functions)  we  consi¬ 
dered.  The  clustering  gets  relatively  stronger  when  p  grows  or  when  the 
side  shape  functions  are  orthogonal ized. 


B.  Multi-p  iterations.  The  multigrid  method  is  widely  used  in  connection 
with  the  finite  difference  method  (see  e.g.  [9]).  The  same  formal  idea  can  be 
used  for  the  p-method  when  exploiting  the  hierarchic  character  of  the  bases 
functions,  see  also  [3]. 

We  will  now  address  the  multi-p  iterative  procedures  also  in  combination 
with  inner  type  iterations.  We  will  address  only  the  bases  functions  of  the 
set  which  were  introduced  in  section  3.  Mostly  we  will  consider  only  the 

case  p  =  8;  the  cases  p  <  8  will  be  mentioned  only  briefly.  First,  let  us 
introduce  some  basic  notions. 

1)  The  hierarchies.  Let  us  describe  the  hierarchies  we  will  later  use. 
The  hierarchy  consists  of  the  sequence  {p^}  of  levels  (degrees)  which  will 
be  iterated  one  at  a  time  in  the  sequel. 


Hj  =  {87654321} 

H2  =  {8421} 

H3  =  {8642} 

H.  =  {87531}. 

4 

By  iteration  on  the  level  p^,  we  mean  iteration  of  all  shape  functions  of 
degree  <  p^  while  the  shape  functions  fo  degree  p  >  p^  are  kept  "frozen. " 

2)  As  in  the  classical  multigrid  we  have  to  deal  with  cycles.  We  will 
use  three  different  cycles. 

\CYCL(Hj.):  It  consists  of  repetition  of  the  hierarchies.  For  example, 
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XCYCLCHj )  =  {876543218765. . .2187. 

/CYCL(K^):  It  consists  of  the  repetition  of  the  hierarchies  in  the 

opposite  order.  For  example  /CYCLCHg)  =  {1248,1248,...}. 

\/CYCL(H^):  It  consists  in  the  repetition  of  the  hierarchies  in  the 
combination  form  \and/.  For  example,  \/CYCL(H_)  =  {8642468642...}. 

On  every  level  p  we  will  perform  s  "inner"  SOR(w)  iterations  with 
the  overrelaxation  parameter  u.  Obviously  w  =  1  coincides  with  the  Gauss 
Seidel  iteration.  The  ordering  is  as  mentioned  in  section  3. 

Table  7.3  shows  the  asymptotic  convergence  rate  p  in  the  H1-norm 
(i.e.,  the  decay  of  the  error  measured  in  the  H^-norm  when  performing  one 
entire  iteration  cycle)  for  p  =  8,  s  =  5,  o>  =  1.8  and  u  =  l. 


Table  7.3 

The  rate  of  convergence  p  for 
various  cycles  and  hierarchies  (p  =  8,  s  *  5). 


\CYCLE 
u  =  1.8 

/CYCLE 
u  =  1.8 

X/CYCLE 
w  =  1.8 

XCYCLE 
u>  =  1.0 

H1 

0.425 

0.537 

0.323 

0.686 

H2 

0.514 

0.507 

0.535 

0.820 

H3 

0.463 

0.  363 

0.497 

0.731 

H4 

0.445 

0.355 

0.449 

0.731 

The  value  u  =  1.8  is  close  to  the  optimal  for  p  =  8  and  all  hierarchies. 
As  ax  illustration  we  show  in  Table  7.4  the  dependence  of  p  on  u>,  for 
\CYCUHj)  and  s  =  5. 
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Table  7.  4. 

The  rate  of  convergence  p  for  \CYCL(H^),  s  =  5,  p  =  8. 


u 

P 

b> 

P 

0 

1.0 

1.6 

0.  446 

0.2 

0.947 

1.7 

0.428 

0.4 

0.891 

1.8 

0.  425 

0.6 

0.829 

1.825 

0.426 

0.8 

0.761 

1.85 

0.435 

1.0 

0.686 

1.875 

0.505 

1.2 

0.603 

1.9 

0.583 

1.4 

0.516 

2.0 

1.00 

Table  7.4  has  shown  the  performances  of  various  cycles  for  u  =  1.0  and  1.8 
and  for  p  =  8.  The  results  for  p  <  8  are  analogous.  As  a  typical  case,  we 
show  the  convergence  rate  for  NCYCLCH^ )  with  s  =  5,  u  =  1.8  and  w  =  1.0. 

Table  7.5. 

The  rate  of  convergence  p  for  NCYCLCH^  with  s  ■  5,  w  *  1,  1.8 


P 

4 

5 

6 

7 

00 

U)  =  1 

0.380 

0.411 

0.  499 

0.516 

0.686 

u  =  1.8 

0.325 

0.323 

0.369 

0.331 

0.425 

So  far  we  presented  the  results  only  for  s  =  5.  In  Table  7.6  we  show  the 
dependence  of  the  rate  p  on  s  in  the  case  XCYCLIH^)  for  p  =  8  and 
w  =  1.8. 
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Table  7.6. 


The  rate  of  convergence  p  for  VCYCLECH^,  u  =  1.8,  p  =  8. 


s 

1 

2 

3 

4 

5 

p 

0.890 

0.748 

0.514 

0.508 

0.  426 

s 

6 

7 

8 

9 

10 

p 

0.379 

0.298 

0.  258 

0.220 

0.  196 

Above  we  have  shown  the  rate  p  for  the  multi-p  iteration.  In  Table  7.7  we 
show  the  rate  p  for  the  SOR  iteration  for  p  =  8. 


Table  7.7. 

The  convergence  rate  p  for  SQR(u),  w  =  1.0,  1.8  and  p  =  8. 


S0R(1.0) 

SOR( 1.8) 

P 

0.966 

0.943 

From  Tables  7.6  and  7.7  we  easily  conclude  that  the  multi-p  iteration  is  more 
effective  (smaller  number  of  operations)  than  the  SOR  iteration  for  p  =  8. 

The  overrelaxation  clearly  improves  the  performance.  The  \CYCL(H^)  is  most 
effective.  Nevertheless  we  conclude  that  the  multi-p  approach  is  less  effec¬ 
tive  than  in  the  classical  multigrid  method  for  the  finite  different  equa¬ 
tions,  at  least  in  the  form  as  we  have  used.  The  section  of  the  shape  func¬ 
tions  could  obviously  also  be  geared  to  the  optimal  performance  of  the  multi-p 
method,  but  the  problem  seems  wide  open. 

We  further  conclude  that  this  conjugate  gradient  method  applied  to  the 
condensed  elements  is  clearly  superior. 

Let  us  now  briefly  compare  the  two  iterative  methods  we  have  addressed. 
The  conjugate  gradient  method  is  very  effective  when  using  the  condensed  ele- 


32 


ments  and  possibly  orthogonal ized  shape  functions.  The  error  decays  (in  the 
energy  norm)  by  one  iteration  by  a  factor  2  or  more  (for  element  with 

p  <  11). 

The  multi-p  method  with  the  hierarchic  shape  functions  does  not  use 

the  condensation  which  was  based  on  the  elimination  procedure  which  is  rela¬ 
tively  expensive  especially  for  high  p.  On  the  other  hand,  it  leads  to  the 
need  of  more  iterations  when  counting  the  inner  iterations  of  the  entire 
cycle.  The  combined  overrelaxation  scheme  behaves  similarly.  Nevertheless, 
counting  all  operations  of  a  cycles  in  comparison  with  the  direct  Gauss-Seidel 
or  overrelaxation  the  multi-p  performance  is  somewhat  better. 

The  conjugate  gradient  method  is  very  good  if  condensation  is  made. 
Without  condensation  the  rate  is  not  good  (e.g. ,  for  p  =  8  we  get  p  = 

0.93).  Among  the  schemes  studied  we  recommend  the  conjugate  gradient  method 


with  condensation  and  with  preconditioning  with  p  =  1  in  the  multielement 
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8.  Conclusions. 

The  selection  of  the  shape  functions  is  of  major  importance  for  the  per¬ 
formance  of  the  solver  based  on  iterative  methods.  Neither  the  theory  nor 
practice  of  the  optimal  selection  of  the  shape  functions  is  available  yet. 

We  have  seen  that  the  condensation  approach  which  has  obvious  advantages 
from  the  point  of  parallel  computations  is  a  very  effective  tool  for  keeping 
the  condition  number  under  the  control  and  is  especially  advantageous  for  the 
conjugate  gradient  method.  The  set  of  the  shape  functions  performs 

better  than  in  the  practical  range  of  p.  Although  we  addressed  only  the 

case  of  one  element  only  a  very  similar  situation  occurs  for  the  mesh  of 
elements  where  p  =  1  preconditioning  is  made.  For  the  analysis  of  the  per¬ 
formance  of  the  conjugate  gradient  method  on  complex  meshes  we  refer  to  [2] 
where  the  parallel  computation  aspects  are  especially  addressed. 

For  the  theoretical  analysis  we  refer  to  [7J. 

We  discussed  only  quadrilateral  elements.  We  expect  that  the  triangular 
elements  will  lead  to  a  similar  performance,  but  extensive  tests  are  not 
available  yet. 

Finally  we  would  like  to  underline  that  we  addressed  only  two-dimensional 
cases.  One  cannot  extrapolate  from  these  results  the  behavior  of  three- 
dimensional  elements. 
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